
/*****************************************************************************************/
/////////////////////////////////// PART FOR STATISTICS ///////////////////////////////////
/*****************************************************************************************/
int     i,t;

// SET STATISTICS TO ZERO //
fracWWtoJ=0.0;fracWWtoUi=0.0;fracUitoWW=0.0;fracUitoJ=0.0;fracJtoU=0.0;fracJtoWW=0.0;fracUntoWW=0.0;fracUntoJ=0.0;necessityUi=0.0;necessityUn=0.0;necessityWW=0.0;*totalassets= 0.0;*totallabor = 0.0;capitalJ=0.0;laborJ=0.0;productionJ=0.0;
fracUi = 0.0;fracUn = 0.0;fracW = 0.0;fracD = 0.0;fracJ = 0.0;frac_employer = 0.0;
expenditure = 0.0;revenubasis = 0.0;corpbase  = 0.0;
welfareUi = 0.0;welfareUn = 0.0;welfareW = 0.0;welfareJ = 0.0;welfaretot = 0.0;C_welfaretot = 0.0;
meanwage = 0.0;avgfirmsize = 0.0;defaultrate = 0.0;costSEA = 0.0;AXtot = 0.0;Atot = 0.0;totalassetsJ = 0.0;leverage= 0.0;
double frac_employer_ability[maxtheta] = {0.0};

for(igrid = 0; igrid < maxgrid; igrid++) {
    distasset[igrid] = 0.0; distassetJ[igrid] = 0.0; distassetW[igrid] = 0.0; distassetUi[igrid] = 0.0; distassetUn[igrid] = 0.0;
    for(tgrid = 0 ; tgrid < maxtheta; tgrid++) {
        fracJpability[tgrid]  = 0.0; fracWWpability[tgrid] = 0.0; fracUpability[tgrid]  = 0.0;
        frac_employer_ability[tgrid] = 0.0;
        necessityUiTHETA[tgrid] = 0.0; necessityUnTHETA[tgrid] = 0.0; necessityWWTHETA[tgrid] = 0.0;
        aggJtoWW[tgrid] = 0.0;aggWWtoJ[tgrid] = 0.0; aggJtoU[tgrid]  = 0.0; aggWWtoU[tgrid] = 0.0; aggUtoWW[tgrid] = 0.0; aggUtoJ[tgrid]  = 0.0;
        ENT[tgrid] = 0.0; SUM[tgrid] = 0.0; AVGFIRM[tgrid] = 0.0; AVGASSET[tgrid] = 0.0; ENTfrac[tgrid]  = 0.0;
        disttotal[inx(igrid, tgrid)] = 0.0;
        for(ygrid = 0; ygrid < maxyprod; ygrid++){distWWearning[inxEARN(tgrid, ygrid)] = 0.0;}
    }
}



/***************************************/
/** START COMPUTE RELEVANT STATISTICS **/
/***************************************/

double incJE, incJNE, fracnegincome, fracnegbusincome, busincJNE, busincJE;
fracnegincome = 0.0; fracnegbusincome = 0.0;

for(igrid=0;igrid<maxgrid;igrid++){
    for(tgrid=0;tgrid<maxtheta;tgrid++){

        /** ENTREPRENEURS' BLOCK **/
        for(zgrid=0;zgrid<maxfirmtype; zgrid++) {


            // COMPUTE TOTAL AX/A, leverage //
            AXtot += (collateralJE[inxE(igrid,tgrid,zgrid)]*distJE[inxE(igrid,tgrid,zgrid)]+collateralJNE[inxE(igrid,tgrid,zgrid)]*distJNE[inxE(igrid,tgrid,zgrid)]);
            Atot += (grid[igrid]*distJE[inxE(igrid,tgrid,zgrid)]+grid[igrid]*distJNE[inxE(igrid,tgrid,zgrid)]);
            if(grid[igrid] > 0){leverage += (sloanJNE[inxE(igrid,tgrid,zgrid)]/grid[igrid])*distJNE[inxE(igrid,tgrid,zgrid)];}

            // POOL OF ENTREPRENEURS //
            ENT[tgrid]      += distJE[inxE(igrid,tgrid,zgrid)]+distJNE[inxE(igrid,tgrid,zgrid)];
            SUM[tgrid]      += distJE[inxE(igrid,tgrid,zgrid)]+distJNE[inxE(igrid,tgrid,zgrid)];
            AVGASSET[tgrid] += mzprob[zgrid][e]*(distJE[inxE(igrid,tgrid,zgrid)]*grid[igrid]+distJNE[inxE(igrid,tgrid,zgrid)]*grid[igrid]);

            // DISTRIBUTION, CAPITAL and WELFARE //
            disttotal[inx(igrid, tgrid)] += distJE[inxE(igrid, tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)];
            fracJ   +=  distJE[inxE(igrid, tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)];
            fracJpability[tgrid] += distJE[inxE(igrid, tgrid, zgrid)]+distJNE[inxE(igrid, tgrid, zgrid)];
            aggJtoWW[tgrid]      += distJEtoWWE[inxE(igrid, tgrid, zgrid)]+distJEtoWWNE[inxE(igrid, tgrid, zgrid)] + distJNEtoWWE[inxE(igrid, tgrid, zgrid)] + distJNEtoWWNE[inxE(igrid, tgrid, zgrid)];
            aggJtoU[tgrid]       += distJEtoUiE[inxE(igrid, tgrid, zgrid)] + distJEtoUnE[inxE(igrid, tgrid, zgrid)] + distJEtoUiNE[inxE(igrid, tgrid, zgrid)] + distJEtoUnNE[inxE(igrid, tgrid, zgrid)] + distJNEtoUiE[inxE(igrid, tgrid, zgrid)] + distJNEtoUnE[inxE(igrid, tgrid, zgrid)] + distJNEtoUiNE[inxE(igrid, tgrid, zgrid)] + distJNEtoUnNE[inxE(igrid, tgrid, zgrid)];

            fracJtoU             += distJEtoUiE[inxE(igrid, tgrid, zgrid)] + distJEtoUnE[inxE(igrid, tgrid, zgrid)] + distJEtoUiNE[inxE(igrid, tgrid, zgrid)] + distJEtoUnNE[inxE(igrid, tgrid, zgrid)] + distJNEtoUiE[inxE(igrid, tgrid, zgrid)] + distJNEtoUnE[inxE(igrid, tgrid, zgrid)] + distJNEtoUiNE[inxE(igrid, tgrid, zgrid)] + distJNEtoUnNE[inxE(igrid, tgrid, zgrid)];
            fracJtoWW            += distJEtoWWE[inxE(igrid, tgrid, zgrid)]+distJEtoWWNE[inxE(igrid, tgrid, zgrid)] + distJNEtoWWE[inxE(igrid, tgrid, zgrid)] + distJNEtoWWNE[inxE(igrid, tgrid, zgrid)];

            distassetJ[igrid]    += distJE[inxE(igrid, tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)];
            totalassetsJ         += (distJE[inxE(igrid, tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)])*grid[igrid];

            welfareJ += distJE[inxE(igrid, tgrid, zgrid)]*valueJE[inxE(igrid, tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)]*valueJNE[inxE(igrid, tgrid, zgrid)]; // Compute the ex-ante welfare of an entrepreneur (before his shock realization).
            C_welfaretot += distJE[inxE(igrid, tgrid, zgrid)]*C_JE[inxE(igrid, tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)]*C_JNE[inxE(igrid, tgrid, zgrid)];
            
            *totalassets += distJE[inxE(igrid, tgrid, zgrid)]*(grid[igrid] - collateralJE[inxE(igrid,tgrid, zgrid)]) + distJNE[inxE(igrid, tgrid, zgrid)]*(grid[igrid]-collateralJNE[inxE(igrid,tgrid, zgrid)]); // compute the K in the corporate sector (- capital used for entrepreneurial sector).

            capitalJ     += distJE[inxE(igrid, tgrid, zgrid)]*collateralJE[inxE(igrid,tgrid, zgrid)] + distJNE[inxE(igrid, tgrid, zgrid)]*(collateralJNE[inxE(igrid,tgrid, zgrid)]); // capital amount used in the entrepreneurial sector


            #if SEA == 1
                
                ENT[tgrid]      += distJEi[inxE(igrid,tgrid,zgrid)]+distJNEi[inxE(igrid,tgrid,zgrid)];
                SUM[tgrid]      += distJEi[inxE(igrid,tgrid,zgrid)]+distJNEi[inxE(igrid,tgrid,zgrid)];
                AVGASSET[tgrid] += distJEi[inxE(igrid,tgrid,zgrid)]*grid[igrid]+distJNEi[inxE(igrid,tgrid,zgrid)]*grid[igrid];
                
                
                disttotal[inx(igrid, tgrid)] += distJEi[inxE(igrid, tgrid, zgrid)] + distJNEi[inxE(igrid, tgrid, zgrid)];
                fracJ                += distJEi[inxE(igrid, tgrid, zgrid)] + distJNEi[inxE(igrid, tgrid, zgrid)];
                fracJi               += distJNEi[inxE(igrid, tgrid, zgrid)] + distJEi[inxE(igrid, tgrid, zgrid)];
                fracJpability[tgrid] += distJEi[inxE(igrid, tgrid, zgrid)]+distJNEi[inxE(igrid, tgrid, zgrid)];
                distassetJ[igrid]    += distJEi[inxE(igrid, tgrid, zgrid)] + distJNEi[inxE(igrid, tgrid, zgrid)];
                welfareJ             += distJEi[inxE(igrid, tgrid, zgrid)]*valueJEi[inxE(igrid, tgrid, zgrid)] + distJNEi[inxE(igrid, tgrid, zgrid)]*valueJNEi[inxE(igrid, tgrid, zgrid)]; // Compute the ex-ante welfare of an entrepreneur (before his shock realization).
                C_welfaretot         += distJEi[inxE(igrid, tgrid, zgrid)]*C_JEi[inxE(igrid, tgrid, zgrid)] + distJNEi[inxE(igrid, tgrid, zgrid)]*C_JNEi[inxE(igrid, tgrid, zgrid)]; // Compute the ex-ante welfare of an entrepreneur (before his shock realization).
                *totalassets         += distJEi[inxE(igrid, tgrid, zgrid)]*(grid[igrid]-collateralJEi[inxE(igrid,tgrid, zgrid)]) + distJNEi[inxE(igrid, tgrid, zgrid)]*(grid[igrid]-(collateralJNEi[inxE(igrid,tgrid, zgrid)]));
                totalassetsJ         += (distJNEi[inxE(igrid, tgrid, zgrid)] + distJEi[inxE(igrid, tgrid, zgrid)])*grid[igrid];
                capitalJ+=distJEi[inxE(igrid, tgrid, zgrid)]*collateralJEi[inxE(igrid,tgrid, zgrid)] + distJNEi[inxE(igrid, tgrid, zgrid)]*(collateralJNEi[inxE(igrid,tgrid, zgrid)]);

            #endif



            // FRACTION OF ENTREPRENEUR WITH NEGATIVE INCOME //
            for(z=0; z<maxfirmtype; z++) {

                // COMPUTE THE PRODUCTION, AVERAGE FIRM SIZE, TAXES, POLICY COST //
                productionJ += (mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)]*prodE(statez[z],tgrid,collateralJE[inxE(igrid,tgrid, zgrid)]) + mzprob[zgrid][z]*distJNE[inxE(igrid, tgrid, zgrid)]*prodE(statez[z],tgrid,collateralJNE[inxE(igrid,tgrid, zgrid)]); // OK
                avgfirmsize += (mzprob[zgrid][z])*(distJE[inxE(igrid, tgrid, zgrid)]*(collateralJE[inxE(igrid,tgrid, zgrid)])+ distJNE[inxE(igrid, tgrid, zgrid)]*(collateralJNE[inxE(igrid,tgrid, zgrid)])); // OK
                AVGFIRM[tgrid] += mzprob[zgrid][z]*(distJE[inxE(igrid,tgrid,zgrid)]*(collateralJE[inxE(igrid,tgrid,zgrid)])+distJNE[inxE(igrid,tgrid,zgrid)]*(collateralJNE[inxE(igrid,tgrid,zgrid)])); // ok

                laborJ+=(mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)]*labor_demand(statez[z],tgrid,collateralJE[inxE(igrid,tgrid, zgrid)]) + mzprob[zgrid][z]*distJNE[inxE(igrid, tgrid, zgrid)]*labor_demand(statez[z],tgrid,collateralJNE[inxE(igrid,tgrid, zgrid)]); // OK
                
                if(labor_demand(statez[z],tgrid,collateralJE[inxE(igrid,tgrid, zgrid)]) > 0){
                    frac_employer += (mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)];
                    frac_employer_ability[tgrid] += (mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)];
                }
                if(labor_demand(statez[z],tgrid,collateralJNE[inxE(igrid,tgrid, zgrid)]) > 0){
                    frac_employer += (mzprob[zgrid][z])*distJNE[inxE(igrid, tgrid, zgrid)];
                    frac_employer_ability[tgrid] +=  (mzprob[zgrid][z])*distJNE[inxE(igrid, tgrid, zgrid)];
                }

                // COMPUTE ENTREPRENEURIAL INCOME //
                if(sloanJNE[inxE(igrid, tgrid, zgrid)] > 0){indicatorAX = 1.0;}else{indicatorAX = 0.0;}
                
                incJE   = (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJE[inxE(igrid, tgrid, zgrid)]) + rstar*(collateralJE[inxE(igrid, tgrid, zgrid)] - grid[igrid]));
                incJNE  = (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJNE[inxE(igrid, tgrid, zgrid)]) - (rateJNE[inxE(igrid, tgrid, zgrid)]-1.0)*(sloanJNE[inxE(igrid, tgrid, zgrid)])) + rstar*(1.0 - indicatorAX)*(grid[igrid]-collateralJNE[inxE(igrid, tgrid, zgrid)]);
                
                busincJE    = (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJE[inxE(igrid, tgrid, zgrid)]));
                busincJNE   = (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJNE[inxE(igrid, tgrid, zgrid)]) - (rateJNE[inxE(igrid, tgrid, zgrid)]-1.0)*(sloanJNE[inxE(igrid, tgrid, zgrid)]));
                
                if(incJE < 0.0){fracnegincome   += distJE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                if(incJNE < 0.0){fracnegincome  += distJNE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                
                if(busincJE < 0.0){fracnegbusincome  += distJE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                if(busincJNE < 0.0){fracnegbusincome += distJNE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                
                
                // COMPUTE TAX BASE + DEFAULT RATE // (in the calibration, corptax == 0). // OK
                returnsbeforetaxJE  = profitE(statez[z],tgrid,collateralJE[inxE(igrid, tgrid, zgrid)]); // OK
                returnsbeforetaxJNE = profitE(statez[z],tgrid,collateralJNE[inxE(igrid, tgrid, zgrid)]) - (rateJNE[inxE(igrid, tgrid, zgrid)])*sloanJNE[inxE(igrid, tgrid, zgrid)]; // OK

                if(returnsbeforetaxJE <= 0.0){basebeforetaxJE = 0.0;}else{basebeforetaxJE = returnsbeforetaxJE;}
                if(returnsbeforetaxJNE <= 0.0){basebeforetaxJNE = 0.0;}else{basebeforetaxJNE = returnsbeforetaxJNE;}


                corpbase += (distJE[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]*basebeforetaxJE + distJNE[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]*basebeforetaxJNE); // OK

                if(defaultJNE[inxEE(igrid, tgrid, zgrid, z)] == 1.0){defaultrate += distJNE[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z];} // OK



                
                #if SEA == 1
                // COMPUTE PRODUCTION  & FIRM SIZE//
                productionJ     += (mzprob[zgrid][z])*distJEi[inxE(igrid, tgrid, zgrid)]*prodE(statez[z],tgrid,collateralJEi[inxE(igrid, tgrid, zgrid)]) + mzprob[zgrid][z]*distJNEi[inxE(igrid, tgrid, zgrid)]*prodE(statez[z],tgrid,collateralJNEi[inxE(igrid, tgrid, zgrid)]); // OK
                avgfirmsize     += (mzprob[zgrid][z])*(distJEi[inxE(igrid, tgrid, zgrid)]*(collateralJEi[inxE(igrid,tgrid, zgrid)])+ distJNEi[inxE(igrid, tgrid, zgrid)]*(collateralJNEi[inxE(igrid,tgrid, zgrid)])); // OK
                AVGFIRM[tgrid]  += mzprob[zgrid][z]*(distJEi[inxE(igrid,tgrid,zgrid)]*(collateralJEi[inxE(igrid,tgrid,zgrid)])+distJNEi[inxE(igrid,tgrid,zgrid)]*(collateralJNEi[inxE(igrid,tgrid,zgrid)])); // ok

                laborJ          +=  (mzprob[zgrid][z])*distJEi[inxE(igrid, tgrid, zgrid)]*labor_demand(statez[z],tgrid,collateralJEi[inxE(igrid,tgrid, zgrid)]) + mzprob[zgrid][z]*distJNEi[inxE(igrid, tgrid, zgrid)]*labor_demand(statez[z],tgrid,collateralJNEi[inxE(igrid,tgrid, zgrid)]); // OK
                
                if(labor_demand(statez[z],tgrid,collateralJEi[inxE(igrid,tgrid, zgrid)]) > 0){
                    frac_employer += (mzprob[zgrid][z])*distJEi[inxE(igrid, tgrid, zgrid)];
                }
                if(labor_demand(statez[z],tgrid,collateralJNEi[inxE(igrid,tgrid, zgrid)]) > 0){
                    frac_employer += (mzprob[zgrid][z])*distJNEi[inxE(igrid, tgrid, zgrid)];
                }
                
                // ENTREPRENEURIAL INCOME //
                if(sloanJNEi[inxE(igrid, tgrid, zgrid)] > 0){indicatorAX = 1.0;}else{indicatorAX = 0.0;}
                
                incJE = (profitE(statez[z],tgrid,collateralJEi[inxE(igrid, tgrid, zgrid)]) + rstar*(collateralJEi[inxE(igrid, tgrid, zgrid)] - grid[igrid]));
                incJNE = (profitE(statez[z],tgrid,collateralJNEi[inxE(igrid, tgrid, zgrid)]) - (rateJNEi[inxE(igrid, tgrid, zgrid)]-1.0)*(collateralJNEi[inxE(igrid, tgrid, zgrid)] - grid[igrid])*indicatorAX) + rstar*(1.0 - indicatorAX)*(grid[igrid]-collateralJNEi[inxE(igrid, tgrid, zgrid)]);
                
                busincJE = (profitE(statez[z],tgrid,collateralJEi[inxE(igrid, tgrid, zgrid)]));
                busincJNE = (profitE(statez[z],tgrid,collateralJNEi[inxE(igrid, tgrid, zgrid)]) - (rateJNEi[inxE(igrid, tgrid, zgrid)]-1.0)*(collateralJNEi[inxE(igrid, tgrid, zgrid)] - grid[igrid])*indicatorAX);
                
                
                if(incJE < 0.0){fracnegincome += distJEi[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                if(incJNE < 0.0){fracnegincome += distJNEi[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                
                if(busincJE < 0.0){fracnegbusincome += distJEi[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}
                if(busincJNE < 0.0){fracnegbusincome += distJNEi[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z];}

                
                // COST OF THE POLICY AND TAXES //
                //returnsbeforetaxJEi  = profitE(statez[z],tgrid,collateralJEi[inxE(igrid, tgrid, zgrid)]);
                //returnsbeforetaxJNEi = profitE(statez[z],tgrid,collateralJNEi[inxE(igrid, tgrid, zgrid)]) - (rateJNEi[inxE(igrid, tgrid, zgrid)])*(collateralJNEi[inxE(igrid, tgrid, zgrid)] - grid[igrid])*indicatorAX;
            
                
                // NOT USED:::
                //if(returnsbeforetaxJEi <= 0.0){basebeforetaxJEi = 0.0;}else{basebeforetaxJEi = returnsbeforetaxJEi;}
                //if(returnsbeforetaxJNEi <= 0.0){basebeforetaxJNEi = 0.0;}else{basebeforetaxJNEi = returnsbeforetaxJNEi;}
                //corpbase += distJEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]*basebeforetaxJEi + distJNEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]*basebeforetaxJNEi; // OK
                // ----
            
                if(defaultJNEi[inxEE(igrid, tgrid, zgrid, z)] == 1.0){
                    defaultrate += distJNEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z];
                    costSEA     += bstarcostfun(0.0, tgrid)*distJNEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]; // OK
                }else{
                    costSEA += bstarcostfun(profitRfun(igrid, z, tgrid, collateralJEi[inxE(igrid, tgrid, zgrid)], rateJNEi[inxE(igrid, tgrid, zgrid)]),tgrid)*distJNEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]; // OK
                }
            
                costSEA += bstarcostfun(profitEhatfun(igrid, z, tgrid, collateralJEi[inxE(igrid, tgrid, zgrid)]), tgrid)*distJEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]; // OK
                
                // COMPUTE DRI COST //
                // not used... corporate taxation.
                //if(returnsbeforetaxJEi <= 0.0){returnsaftertaxJEi = returnsbeforetaxJEi;}else{returnsaftertaxJEi = returnsbeforetaxJEi*(1 - corptax);} // OK
            
                
            
                // for JNE compute ///
                // not used... corporate taxation.
//                if(corptax > 0.0){printf("corptax is not coded at the moment.");getchar();}
//                if(defaultJNEi[inxEE(igrid, tgrid, zgrid, z)] == 0){
//                    if(returnsbeforetaxJNEi <= 0.0){returnsaftertaxJNEi = returnsbeforetaxJNEi;}else{returnsaftertaxJNEi = returnsbeforetaxJNEi*(1 - corptax);} // OK
//                    costSEA += bstarcostfun(returnsaftertaxJNEi, tgrid)*distJNEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]; // OK
//                    
//                } else {
//                    costSEA += bstarcostfun(0.0, tgrid)*distJNEi[inxE(igrid, tgrid, zgrid)]*mzprob[zgrid][z]; // OK
//                }
                #endif
            

            } // end z
        } // end zgrid




        /** BLOCK WORKER **/
        for(ygrid = 0; ygrid < maxyprod; ygrid++) {

            // DISTRIBUTION //
            SUM[tgrid]                   += distWWNE[inxW(igrid,tgrid, ygrid)] + distWWE[inxW(igrid,tgrid, ygrid)];
            disttotal[inx(igrid, tgrid)] += distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)];
            fracWWpability[tgrid]        += distWWNE[inxW(igrid, tgrid, ygrid)] + distWWE[inxW(igrid, tgrid, ygrid)];
            aggWWtoJ[tgrid]              += distWWEtoJE[inxW(igrid, tgrid, ygrid)] + distWWEtoJNE[inxW(igrid, tgrid, ygrid)] + distWWNEtoJNE[inxW(igrid, tgrid, ygrid)];
            aggWWtoU[tgrid]              += distWWEtoUiNE[inxW(igrid, tgrid, ygrid)] + distWWEtoUiE[inxW(igrid, tgrid, ygrid)] + distWWNEtoUiNE[inxW(igrid, tgrid, ygrid)];
            fracWWtoJ                    += distWWEtoJE[inxW(igrid, tgrid, ygrid)] + distWWEtoJNE[inxW(igrid, tgrid, ygrid)] + distWWNEtoJNE[inxW(igrid, tgrid, ygrid)];
            fracWWtoUi                   += distWWEtoUiNE[inxW(igrid, tgrid, ygrid)] + distWWEtoUiE[inxW(igrid, tgrid, ygrid)] + distWWNEtoUiNE[inxW(igrid, tgrid, ygrid)];
            distassetW[igrid]            += distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)];
            distWWearning[inxEARN(tgrid, ygrid)] += distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)];
            // K //
            *totalassets            += (distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)])*grid[igrid];
            // L //
            *totallabor             += (distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)])*prod[tgrid]*yprod[ygrid];
            // FRAC WW //
            fracW                   += distWWNE[inxW(igrid, tgrid, ygrid)] + distWWE[inxW(igrid, tgrid, ygrid)]; // OK
            // WELFARE //
            welfareW                += distWWE[inxW(igrid, tgrid, ygrid)]*valueWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)]*valueWWNE[inxW(igrid, tgrid, ygrid)]; // OK
            C_welfaretot            += distWWE[inxW(igrid, tgrid, ygrid)]*C_WWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)]*C_WWNE[inxW(igrid, tgrid, ygrid)]; // OK
            // NECESSITY SHARE //
            necessityWW             += necessityWWE[inxW(igrid, tgrid, ygrid)] + necessityWWNE[inxW(igrid, tgrid, ygrid)]; // OK
            necessityWWTHETA[tgrid] += necessityWWE[inxW(igrid, tgrid, ygrid)] + necessityWWNE[inxW(igrid, tgrid, ygrid)]; // OK
            // FOR TAXES // OK
            revenubasis             += (distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)])*wstar*prod[tgrid]*yprod[ygrid];
            // MEAN WAGE //
            meanwage                += (distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)])*yprod[ygrid]*prod[tgrid]*wstar*(1.0 - taxrate);

            // EARNING DISTRIBUTION //
            earningdist[inxW(igrid, tgrid, ygrid)][0] = prod[tgrid]*yprod[ygrid]*wstar*(1.0 - taxrate) + 0.00000001*igrid;
            earningdist[inxW(igrid, tgrid, ygrid)][1] = distWWNE[inxW(igrid, tgrid, ygrid)] + distWWE[inxW(igrid, tgrid, ygrid)];
            earningdist[inxW(igrid, tgrid, ygrid)][2] = distWWEtoJE[inxW(igrid, tgrid, ygrid)] + distWWEtoJNE[inxW(igrid, tgrid, ygrid)] + distWWNEtoJNE[inxW(igrid, tgrid, ygrid)];
            earningdist[inxW(igrid, tgrid, ygrid)][3] = distWWEtoUiNE[inxW(igrid, tgrid, ygrid)] + distWWEtoUiE[inxW(igrid, tgrid, ygrid)] + distWWNEtoUiNE[inxW(igrid, tgrid, ygrid)];

        }


        /** UNEMPLOYED BLOCK **/
        // DISTRIBUTION //
        SUM[tgrid]                  += distUiNE[inx(igrid,tgrid)] + distUiE[inx(igrid,tgrid)] + distUnNE[inx(igrid,tgrid)] + distUnE[inx(igrid,tgrid)];
        *totalassets                += (distUiE[inx(igrid,tgrid)]* + distUnE[inx(igrid,tgrid)] + distUiNE[inx(igrid,tgrid)] + distUnNE[inx(igrid,tgrid)])*grid[igrid];
        disttotal[inx(igrid,tgrid)] += distUiE[inx(igrid, tgrid)] + distUnE[inx(igrid, tgrid)] + distUiNE[inx(igrid, tgrid)] + distUnNE[inx(igrid, tgrid)];
        distassetUi[igrid]          += distUiE[inx(igrid, tgrid)] + distUiNE[inx(igrid, tgrid)];
        distassetUn[igrid]          += distUnE[inx(igrid, tgrid)] + distUnNE[inx(igrid, tgrid)];

        aggUtoJ[tgrid] += distUnEtoJE[inx(igrid, tgrid)] + distUnEtoJNE[inx(igrid, tgrid)] + distUiEtoJE[inx(igrid, tgrid)] + distUiEtoJNE[inx(igrid, tgrid)] + distUnNEtoJNE[inx(igrid, tgrid)] + distUiNEtoJNE[inx(igrid, tgrid)]; // ok
        aggUtoWW[tgrid] += distUnEtoWWNE[inx(igrid, tgrid)] + distUnEtoWWE[inx(igrid, tgrid)] + distUiEtoWWE[inx(igrid, tgrid)] + distUiEtoWWNE[inx(igrid, tgrid)] + distUnNEtoWWNE[inx(igrid, tgrid)] + distUiNEtoWWNE[inx(igrid, tgrid)]; // ok
        fracUitoJ  += distUiEtoJE[inx(igrid, tgrid)] + distUiEtoJNE[inx(igrid, tgrid)] + distUiNEtoJNE[inx(igrid, tgrid)]; // ok
        fracUitoWW += distUiEtoWWE[inx(igrid, tgrid)] + distUiEtoWWNE[inx(igrid, tgrid)] + distUiNEtoWWNE[inx(igrid, tgrid)]; // ok
        fracUntoJ  += distUnEtoJE[inx(igrid, tgrid)] + distUnEtoJNE[inx(igrid, tgrid)] + distUnNEtoJNE[inx(igrid, tgrid)]; // ok
        fracUntoWW += distUnEtoWWNE[inx(igrid, tgrid)] + distUnEtoWWE[inx(igrid, tgrid)] + distUnNEtoWWNE[inx(igrid, tgrid)]; // ok
        fracUpability[tgrid] += distUiNE[inx(igrid, tgrid)] + distUiE[inx(igrid, tgrid)] +distUnNE[inx(igrid, tgrid)] + distUnE[inx(igrid, tgrid)];
        distasset[igrid]     += disttotal[inx(igrid, tgrid)];

        // FRACTION UNEMPLOYED //
        fracUi += distUiE[inx(igrid, tgrid)] + distUiNE[inx(igrid, tgrid)];  // insured
        fracUn += distUnE[inx(igrid, tgrid)] + distUnNE[inx(igrid, tgrid)];  // uninsured

        // WELFARE //
        welfareUi += distUiE[inx(igrid, tgrid)]*valueUiE[inx(igrid,tgrid)] + distUiNE[inx(igrid, tgrid)]*valueUiNE[inx(igrid,tgrid)]; // ok
        welfareUn += distUnE[inx(igrid, tgrid)]*valueUnE[inx(igrid,tgrid)] + distUnNE[inx(igrid, tgrid)]*valueUnNE[inx(igrid,tgrid)]; // ok
        C_welfaretot += distUiE[inx(igrid, tgrid)]*C_UiE[inx(igrid,tgrid)] + distUiNE[inx(igrid, tgrid)]*C_UiNE[inx(igrid,tgrid)]; // ok
        C_welfaretot += distUnE[inx(igrid, tgrid)]*C_UnE[inx(igrid,tgrid)] + distUnNE[inx(igrid, tgrid)]*C_UnNE[inx(igrid,tgrid)]; // ok

        // Necessity share // ok
        necessityUi             +=  necessityUiE[inx(igrid, tgrid)] + necessityUiNE[inx(igrid, tgrid)];
        necessityUn             +=  necessityUnE[inx(igrid, tgrid)] + necessityUnNE[inx(igrid, tgrid)];
        necessityUiTHETA[tgrid] +=  necessityUiE[inx(igrid, tgrid)] + necessityUiNE[inx(igrid, tgrid)];
        necessityUnTHETA[tgrid] +=  necessityUnE[inx(igrid, tgrid)] + necessityUnNE[inx(igrid, tgrid)];

        // UI cost //
        expenditure += (distUiE[inx(igrid,tgrid)] + distUiNE[inx(igrid,tgrid)])*(wstar*rhostar*prod[tgrid]*(1.0-taxrate));
        //expenditure += (distUnE[inx(igrid,tgrid)] + distUnNE[inx(igrid,tgrid)] + distUn[inx(igrid,tgrid)])*minpar; // minimum subsitance level.

    }
}


// NORMALIZE THE VALUES //
for(tgrid = 0; tgrid < maxtheta; tgrid++){
    ENTfrac[tgrid]  = ENT[tgrid]/SUM[tgrid];        // ok
    AVGFIRM[tgrid]  = AVGFIRM[tgrid]/ENT[tgrid];    // ok
    AVGASSET[tgrid] = AVGASSET[tgrid]/ENT[tgrid];   // ok
}

ratioAXA = AXtot/Atot;


frac_employer_ability[0] = frac_employer_ability[0]/ENT[0];
frac_employer_ability[1] = frac_employer_ability[1]/ENT[1];
frac_employer_ability[2] = frac_employer_ability[2]/ENT[2];



/** WAGE **/
meanwage = meanwage/fracW;

/** WELFARE **/
welfaretot  = welfareJ + welfareUn + welfareUi + welfareW;      // ex-ante welfare, before shock on entrepreneurs.
welfareJ    = welfareJ/fracJ;
welfareW    = welfareW/fracW;
welfareUi   = welfareUi/fracUi;
welfareUn   = welfareUn/fracUn;

/** AVERAGE FIRM SIZE **/
avgfirmsize = avgfirmsize/fracJ;

// FOR PM METHOD //
#if indexPM == 1
    *fracJJ         = fracJ;
    *fracUU         = fracUn + fracUi;
    *meanwagelike   = meanwage;
#endif


/** COMPUTE AGGREGATES **/
//*totalassets = *totalassets - defaultcost; // compte the net amount of asset in the economy, taking into account losses from default.
totalcapital    = *totalassets+capitalJ;
corporate_labor = *totallabor-laborJ;
corproduction   = TFPpar*pow((*totalassets),(alphapar))*pow((corporate_labor), (1.0-alphapar))*4.0;  // this is an annual production
totalproduction = (corproduction+productionJ*4); // transform into quarter production.
KY = totalcapital/totalproduction; // OK




/**********************************************/
//      GINI AND WEALTH / INCOME DISTRI       //
/**********************************************/


// COMPUTE DIFFERENT TOP BOTTOM OF THE DISTRIBUTION
GINIval = GINI(grid, distasset, maxgrid);
top001p = 1 - toppercent(grid, distasset, maxgrid, 0.999);
top01p = 1 - toppercent(grid, distasset, maxgrid, 0.99);
top05p = 1 - toppercent(grid, distasset, maxgrid, 0.95);
top10p = 1 - toppercent(grid, distasset, maxgrid, 0.90);
top20p = 1 - toppercent(grid, distasset, maxgrid, 0.80);
top40p = 1 - toppercent(grid, distasset, maxgrid, 0.60);
bot20p = toppercent(grid, distasset, maxgrid, 0.20);
bot40p = toppercent(grid, distasset, maxgrid, 0.40);
bot60p = toppercent(grid, distasset, maxgrid, 0.60);
bot80p = toppercent(grid, distasset, maxgrid, 0.80);



// COMPUTE MEDIAN WORTH
if (fracJ>0.0){medianworth(distassetJ,fracJ, maxgrid, &mworthJ);}else{mworthJ=0.0;}
if (fracW>0.0){medianworth(distassetW,fracW, maxgrid, &mworthW);}else{mworthW=0.0;}
if (fracUi>0.0){medianworth(distassetUi,fracUi, maxgrid, &mworthUi);}else{mworthUi=0.0;}
if (fracUn>0.0){medianworth(distassetUn,fracUn, maxgrid, &mworthUn);}else{mworthUn=0.0;}
medianworth(distasset,1.0, maxgrid, &mworthall);



// COMPUTE EARNINGS DISTRIBUTION (after tax) (gini around 0.36)
double EARNINGS[ifulldimEARN];
for(tgrid=0; tgrid < maxtheta; tgrid++){
    for(ygrid=0; ygrid < maxyprod; ygrid++) {
        EARNINGS[inxEARN(tgrid, ygrid)] = prod[tgrid]*yprod[ygrid]*wstar*(1.0 - taxrate);
    }
}
GINIearn = GINI(EARNINGS, distWWearning, ifulldimEARN);


// COMPUTE INCOME DISTRIBUTION
double varbusinc = 0.0, meanbusinc= 0.0,earningpos = 0.0,stdentinc= 0.0,  varworkinc= 0.0, meanworkinc= 0.0,stdworkinc= 0.0,nb_work= 0.0,fracJ_pos = 0.0;
for(igrid=0;igrid<maxgrid;igrid++){
    for(tgrid=0;tgrid<maxtheta;tgrid++){
        for(zgrid=0;zgrid<maxfirmtype; zgrid++) {
            for(z=0; z<maxfirmtype; z++) {
                earningpos = ((1.0 - corptax)*(profitE(statez[z],tgrid,collateralJE[inxE(igrid, tgrid, zgrid)])));
                if(earningpos > 0.00001){
                    varbusinc  += log(earningpos)*log(earningpos)*(mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)];
                    meanbusinc += log(earningpos)*(mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)];
                    fracJ_pos += (mzprob[zgrid][z])*distJE[inxE(igrid, tgrid, zgrid)];
                }
                earningpos = (((1.0 - corptax)*(profitE(statez[z],tgrid,collateralJNE[inxE(igrid, tgrid, zgrid)]) - (rateJNE[inxE(igrid, tgrid, zgrid)]-1.0)*(sloanJNE[inxE(igrid, tgrid, zgrid)]))));
                if(earningpos > 0.00001){
                    varbusinc  += log(earningpos)*log(earningpos)*(mzprob[zgrid][z])*distJNE[inxE(igrid, tgrid, zgrid)];
                    meanbusinc += log(earningpos)*(mzprob[zgrid][z])*distJNE[inxE(igrid, tgrid, zgrid)];
                    fracJ_pos  += (mzprob[zgrid][z])*distJNE[inxE(igrid, tgrid, zgrid)];
                }
            }
            
        }
        
        for(ygrid = 0; ygrid < maxyprod; ygrid++) {
            varworkinc += log(prod[tgrid]*yprod[ygrid]*wstar*(1.0 - taxrate))*log(prod[tgrid]*yprod[ygrid]*wstar*(1.0 - taxrate))*(distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)]);
            meanworkinc += log(prod[tgrid]*yprod[ygrid]*wstar*(1.0 - taxrate))*(distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)]);
            nb_work += (distWWE[inxW(igrid, tgrid, ygrid)] + distWWNE[inxW(igrid, tgrid, ygrid)]);
        }
    }
}
stdentinc  = pow(varbusinc/fracJ_pos - pow(meanbusinc/fracJ_pos,2.0),0.5);
stdworkinc = pow(varworkinc/nb_work - pow(meanworkinc/nb_work,2.0),0.5);

//printf("%f %f %f %f \n",varbusinc,fracJ,meanbusinc,stdentinc);
//printf("%f %f %f %f \n",varworkinc,nb_work,meanworkinc,stdworkinc); //getchar();


/** COMPUTE THREE BINS OF EARNINGS (THREE BINS) **/

// 1. SORT EARNING DISTRIBUTION BY EARNINGS
qsort(earningdist, ifulldimW, 4*sizeof(double), comparefun2);

// 2. SAVE VALUES OF INDEX FOR WHICH DISTRI <= 0.33 FOR EX
double stop1, stop2, masscumul;
int indexseuil[2];
int l;
masscumul = 0.0; stop1 = 0.0; stop2 = 0.0;
for(l=0;l<ifulldimW;l++){
    masscumul += earningdist[l][1]/fracW;
    if(masscumul >= 0.333  && stop1 == 0.0){indexseuil[0] = l;stop1 = 1.0;}
    if(masscumul >= 0.667 && stop2 == 0.0){indexseuil[1] = l;stop2 = 1.0;}
}


// 3. COMPUTE THE TRANSITION WW to E FOR SUCH INDEXES, how ? NEED to compare all indexes.
double transiWWtoE[3];
transiWWtoE[0] = 0.0;
transiWWtoE[1] = 0.0;
transiWWtoE[2] = 0.0;
for(l=0;l<ifulldimW;l++) {
    if(l <= indexseuil[0]){transiWWtoE[0] += (earningdist[l][2])/(0.333*fracW);}
    if(l > indexseuil[0] && l <= indexseuil[1]){transiWWtoE[1] += (earningdist[l][2])/(0.334*fracW);}
    if(l > indexseuil[1]){transiWWtoE[2] += (earningdist[l][2])/(0.333*fracW);}
}

// 3. COMPUTE THE TRANSITION WW to E FOR SUCH INDEXES, how ? NEED to compare all indexes.
double transiWWtoU[3];
transiWWtoU[0] = 0.0;
transiWWtoU[1] = 0.0;
transiWWtoU[2] = 0.0;
for(l=0;l<ifulldimW;l++) {
    if(l <= indexseuil[0]){transiWWtoU[0] += (earningdist[l][3])/(0.333*fracW);}
    if(l > indexseuil[0] && l <= indexseuil[1]){transiWWtoU[1] += (earningdist[l][3])/(0.334*fracW);}
    if(l > indexseuil[1]){transiWWtoU[2] += (earningdist[l][3])/(0.333*fracW);}
}






/**********************************************/
//                  GOVERNMENT                //
/**********************************************/

#if NOPOLICY == 1
    *optimaltax = (expenditure + GratioY*(corproduction + productionJ*4)/4 - corpbase*corptax)/revenubasis;

    printf("TAX : Optitax: %f Exp: %f  Revbase: %f CorpRev: %f Corpbase: %f Meanwage: %f ProductionJ: %f\n", *optimaltax, expenditure + GratioY*corproduction, revenubasis, corpbase*corptax, corpbase, meanwage, productionJ*4);
#endif

#if SEA == 1
    #if earningSEAtax == 1 // use earning tax rate
            *optimaltaxSEA = lumpsumSEA;
            *optimaltax = (expenditure + costSEA + GratioY*(corproduction + productionJ*4)/4 - corpbase*corptax)/revenubasis;
    #endif
    #if earningSEAtax == 0 // use lump-sum
            *optimaltax = (expenditure + GratioY*(corproduction + productionJ*4)/4 - corpbase*corptax)/revenubasis;
            *optimaltaxSEA = (costSEA)/(fracJ + fracW + fracUi); // if lumpsum tax
    #endif
    printf("TAX : Optitax: %f OptiSEA: %f Exp: %f  Revbase: %f CorpRev: %f Corpbase: %f CostSEA: %f Meanwage: %f ProductionJ: %f fracJi: %f\n", *optimaltax, *optimaltaxSEA, expenditure + GratioY*corproduction, revenubasis, corpbase*corptax, corpbase, costSEA, meanwage, productionJ*4, fracJi);
#endif

#if indexPM == 1
    *defaultratetot = defaultrate/fracJ;
    *prodtotal      = totalproduction;
    *probleaveU     = (fracUntoJ + fracUitoJ)/(fracUi+fracUn);
    *newENTq        = (fracUitoJ + fracUntoJ + fracWWtoJ)/(fracJ+fracW+fracUi+fracUn);
#endif





/**********************************************/
//             PRINT ALL STATS                //
/**********************************************/

tempfileoutfile=fopen(tempfileout, "a");
setbuf ( tempfileoutfile , NULL );

fprintf(tempfileoutfile,"-------EQUILIBRIUM PRICES AND AGGREGATE--------\n");
fprintf(tempfileoutfile,"r %20.15f \n", rstar);
fprintf(tempfileoutfile,"w %20.15f \n", wstar);
fprintf(tempfileoutfile,"mw %20.15f\n", meanwage);
fprintf(tempfileoutfile,"u %20.15f \n", fracUi + fracUn);
#if SEA == 1
            fprintf(tempfileoutfile,"taxSEA %20.15f \n", lumpsumSEA);
#endif

fprintf(tempfileoutfile,"-------MOMENTS GENERATED & STATISTICS--------\n");
fprintf(tempfileoutfile,"Ratio KY annual [12.16] %20.15f\n",KY*4);
fprintf(tempfileoutfile,"Fraction of E>nonE [6.7] %20.15f\n",100*(fracJtoU + fracJtoWW)/fracJ);
fprintf(tempfileoutfile,"Fraction of E [8.2] %20.15f\n",100*fracJ);
fprintf(tempfileoutfile,"Ratio mworthJ/mworthW [8] %20.15f\n",(mworthJ/mworthW));
fprintf(tempfileoutfile,"Wealth Gini [0.8] %20.15f\n",GINIval);
fprintf(tempfileoutfile,"Bankruptcy rate [0.75] %20.15f\n",100*(defaultrate/fracJ));
fprintf(tempfileoutfile,"U>E [2.2] %20.15f\n",100*(fracUntoJ + fracUitoJ)/(fracUi+fracUn));
fprintf(tempfileoutfile,"(UL: %20.15f, US: %20.15f) Unemploy. rate [6.4] %20.15f\n",fracUn, fracUi, 100*(fracUi + fracUn));
fprintf(tempfileoutfile,"IUR (Ui/W) [2.5] = %f\n", (fracUi/(fracW))*100);
fprintf(tempfileoutfile,"W>E 1 [1.07] %20.15f\n",transiWWtoE[0]/((fracWWtoJ)/fracW));
fprintf(tempfileoutfile,"W>E 2 [0.89] %20.15f\n",transiWWtoE[1]/((fracWWtoJ)/fracW));
fprintf(tempfileoutfile,"W>E 3 [1.07] %20.15f\n",transiWWtoE[2]/((fracWWtoJ)/fracW));
fprintf(tempfileoutfile,"W>U 1 [1.07] %20.15f\n",transiWWtoU[0]/((fracWWtoUi)/fracW));
fprintf(tempfileoutfile,"W>U 2 [0.89] %20.15f\n",transiWWtoU[1]/((fracWWtoUi)/fracW));
fprintf(tempfileoutfile,"W>U 3 [1.07] %20.15f\n",transiWWtoU[2]/((fracWWtoUi)/fracW));
fprintf(tempfileoutfile,"Zero wealth fraction [7- 13] %20.15f\n",100*(distasset[0]));


fprintf(tempfileoutfile,"------- OTHER MOMENTS GENERATED & STATISTICS--------\n");
if ((fracJ+fracW+fracUi+fracUn)>0.0){fprintf(tempfileoutfile,"Fraction of nonE>E (rlv to whole pop) [0.39] %20.15f\n",((fracUitoJ + fracUntoJ + fracWWtoJ)/(fracJ+fracW+fracUi+fracUn)));}else{fprintf(tempfileoutfile,"Fraction of nonE>E (rlv to whole pop) %20s\n","nan");} // 5th moment
if ((fracJ+fracW+fracUi+fracUn)>0.0){fprintf(tempfileoutfile,"Fraction of U>E (rlv to new entrant) [20] %20.15f\n",(fracUitoJ + fracUntoJ)/(fracUitoJ + fracUntoJ + fracWWtoJ));}else{fprintf(tempfileoutfile,"Fraction of U>E (rlv to new entrant) [20] %20s\n","nan");} // 6th moment
fprintf(tempfileoutfile,"Worker exit rate (rlv to workers) [3-6.5] %20.15f\n",(fracWWtoJ + fracWWtoUi)/(fracW));
fprintf(tempfileoutfile,"Nshare (US) %20.15f, \t Nshare(UL): %20.15f, \t Nshare(WW): %20.15f, \t Nshare [0.1]: %20.15f\n",(necessityUi)/(fracUitoJ + fracUntoJ + fracWWtoJ), necessityUn/(fracUitoJ + fracUntoJ + fracWWtoJ), necessityWW/(fracUitoJ + fracUntoJ + fracWWtoJ), (necessityUn + necessityUi + necessityWW)/(fracUitoJ + fracUntoJ + fracWWtoJ)); // 11th moment
fprintf(tempfileoutfile,"Earning Gini [0.36] %20.15f\n",GINIearn);                      // 7bisth moment
fprintf(tempfileoutfile,"Top 1 [30] %20.15f\n",top01p);                               // 8th moment
if ((fracUi+fracUn)>0.0){fprintf(tempfileoutfile,"Fraction U>nonU (rlv to unemployed) [?] %20.15f\n",100*(fracUitoJ + fracUitoWW + fracUntoWW + fracUntoJ)/(fracUi+fracUn));}else{fprintf(tempfileoutfile,"Fraction U>nonU out %20s\n","nan");}                  // 10th moment
fprintf(tempfileoutfile,"\n");


fprintf(tempfileoutfile,"-------TRANSITION Worker to entrepreneur / ABILITY--------\n");
fprintf(tempfileoutfile,"very low\t middle \t very high\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",aggWWtoJ[tgrid]/fracWWpability[tgrid]);
}
fprintf(tempfileoutfile,"\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",transiWWtoE[tgrid]);
}
fprintf(tempfileoutfile,"\n");


fprintf(tempfileoutfile,"-------TRANSITION Worker to Unemployed / ABILITY--------\n");
fprintf(tempfileoutfile,"very low \t middle \t very high\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",aggWWtoU[tgrid]/fracWWpability[tgrid]);
}
fprintf(tempfileoutfile,"\n");

fprintf(tempfileoutfile,"-------TRANSITION Entrepreneur to Worker / ABILITY--------\n");
fprintf(tempfileoutfile,"very low \t middle \t very high\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",aggJtoWW[tgrid]/fracJpability[tgrid]);
}
fprintf(tempfileoutfile,"\n");

fprintf(tempfileoutfile,"-------TRANSITION Entrepreneur to Unemployed / ABILITY--------\n");
fprintf(tempfileoutfile,"very low \t middle \t very high\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",aggJtoU[tgrid]/fracJpability[tgrid]);
}
fprintf(tempfileoutfile,"\n");

fprintf(tempfileoutfile,"-------TRANSITION Unemployed to Worker / ABILITY--------\n");
fprintf(tempfileoutfile,"very low \t middle \t very high\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",aggUtoWW[tgrid]/fracUpability[tgrid]);
}
fprintf(tempfileoutfile,"\n");

fprintf(tempfileoutfile,"-------TRANSITION Unemployed to Entrepreneur / ABILITY--------\n");
fprintf(tempfileoutfile,"very low \t middle \t very high\n");
for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
fprintf(tempfileoutfile,"%20.15f\t",aggUtoJ[tgrid]/fracUpability[tgrid]);
}
fprintf(tempfileoutfile,"\n");


fprintf(tempfileoutfile,"--------- TRANSITION MATRIX BTW OCCUPATIONS ------\n");
fprintf(tempfileoutfile, "  \t  U     \t W     \t E \n");
fprintf(tempfileoutfile, "U \t %3.7f \t %3.7f \t %3.7f \n", (1- (fracUitoJ + fracUitoWW + fracUntoWW + fracUntoJ)/(fracUi+fracUn)), (fracUitoWW + fracUntoWW)/(fracUi+fracUn), (fracUntoJ + fracUitoJ)/(fracUi+fracUn));
fprintf(tempfileoutfile, "W \t %3.7f \t %3.7f \t %3.7f \n", (fracWWtoUi)/fracW, (1 - (fracWWtoJ + fracWWtoUi)/fracW), (fracWWtoJ)/fracW);
fprintf(tempfileoutfile, "E \t %3.7f \t %3.7f \t %3.7f \n", (fracJtoU)/fracJ, fracJtoWW/fracJ, (1-(fracJtoWW + fracJtoU)/fracJ));
fprintf(tempfileoutfile,"---------------\n");


fprintf(tempfileoutfile,"--------- WEALTH DISTRIBUTION STATISTICS ------\n");
fprintf(tempfileoutfile,"Top 0.1 [22] %20.15f\n",top001p);
fprintf(tempfileoutfile,"Top 1 [41.8] %20.15f\n",top01p); // 8th moment
fprintf(tempfileoutfile,"Top 5 [?] %20.15f\n",top05p);
fprintf(tempfileoutfile,"Top 10 [77.2] %20.15f\n",top10p);
fprintf(tempfileoutfile,"Top 20 [?] %20.15f\n",top20p);
fprintf(tempfileoutfile,"Top 40 [?] %20.15f\n",top40p);
fprintf(tempfileoutfile,"Bot 20 [?] %20.15f\n",bot20p);
fprintf(tempfileoutfile,"Bot 40 [?] %20.15f\n",bot40p);
fprintf(tempfileoutfile,"Bot 60 [?] %20.15f\n",bot60p);
fprintf(tempfileoutfile,"Bot 80 [?] %20.15f\n",bot80p);
fprintf(tempfileoutfile,"mworthJ %20.15f\n",mworthJ);
fprintf(tempfileoutfile,"mworthW %20.15f\n",mworthW);
fprintf(tempfileoutfile,"mworthUi %20.15f\n",mworthUi);
fprintf(tempfileoutfile,"mworthUn %20.15f\n",mworthUn);
fprintf(tempfileoutfile,"avgfirm: %20.15f\n",avgfirmsize);
fprintf(tempfileoutfile,"---------------\n");


fprintf(tempfileoutfile,"--------- PRODUCTION & CAPITAL -----\n");
fprintf(tempfileoutfile,"capitalJ %20.15f\n",capitalJ);
fprintf(tempfileoutfile,"productionJ (annual) %20.15f\n",productionJ*4);
fprintf(tempfileoutfile,"Totalcapital %20.15f\n",totalcapital);
fprintf(tempfileoutfile,"Totallabor %20.15f\n",*totallabor);
fprintf(tempfileoutfile,"Corproduction (annual) %20.15f\n",corproduction);
fprintf(tempfileoutfile,"Corcapital %20.15f\n",*totalassets);
fprintf(tempfileoutfile,"Totalproduction (annual) %20.15f\n",totalproduction);
fprintf(tempfileoutfile,"capitalJ/Totalcapital %20.15f\n",(capitalJ/totalcapital));
fprintf(tempfileoutfile,"---------------\n");


fprintf(tempfileoutfile,"--------- TAXATION -----\n");
fprintf(tempfileoutfile,"Corporate revenu %20.15f\n",corpbase*corptax);
fprintf(tempfileoutfile,"Expenditure %20.15f\n",expenditure);
fprintf(tempfileoutfile,"Revenue Basis %20.15f\n",revenubasis);
fprintf(tempfileoutfile,"Optimal tax %20.15f\n",*optimaltax);
#if SEA == 1
    fprintf(tempfileoutfile,"Optimal tax SEA %20.15f\n",*optimaltaxSEA);
    fprintf(tempfileoutfile,"COST SEA: %20.15f %20.15f\n",costSEA/(productionJ*4),costSEA/totalproduction);
    fprintf(tempfileoutfile,"COST SEA: %20.15f %20.15f\n",costsubsidy/(productionJ*4),costsubsidy/totalproduction);
#endif

fprintf(tempfileoutfile,"-------MASS--------\n");
fprintf(tempfileoutfile,"WW %20.15f \n", fracW);
fprintf(tempfileoutfile,"U %20.15f \n", fracUn + fracUi);
fprintf(tempfileoutfile,"VE %20.15f \n", fracJ);

fprintf(tempfileoutfile,"-------WELFARE--------\n");
fprintf(tempfileoutfile,"Un %20.15f  Ui %20.15f  J %20.15f  WW %20.15f  Wtot%20.15f   Ctot%20.15f\n", welfareUn, welfareUi, welfareJ, welfareW, welfaretot, C_welfaretot);

fprintf(tempfileoutfile,"-------POOL OF ENT--------\n");
for(tgrid = 0; tgrid < maxtheta; tgrid++) {
fprintf(tempfileoutfile,"THETA: %d,  %20.15f, %20.15f, %20.15f, %20.15f  \n", tgrid, ENTfrac[tgrid], SUM[tgrid], AVGFIRM[tgrid], AVGASSET[tgrid]);
}

fclose(tempfileoutfile);



////////////////////////
// PRINT MOMENTS UiED //
////////////////////////

distoutfile=fopen(moment, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","KY", "E", "U", "TOP1", "GINI", "netVEWW", "ZERO", "NECESSITY","EEXIT", "UEXIT", "UtoE","WWEXIT");
fprintf(distoutfile,"%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",KY, fracJ, fracUi + fracUn, top01p, GINIval, (mworthJ/mworthW), distasset[0], (necessityUn + necessityUi)/(fracUntoJ + fracUitoJ + fracWWtoJ), (fracJtoU + fracJtoWW)/fracJ, (fracUitoWW + fracUitoJ + fracUntoWW + fracUntoJ)/(fracUn + fracUi), (fracUitoJ + fracUntoJ)/(fracUi + fracUn), (fracWWtoJ + fracWWtoUi)/(fracW));
fclose(distoutfile);



/**** PRINT RESULTS ON THE PROGRAM ****/
#if indexPM == 0
    printf("-------MOMENTS GENERATED --------\n");
    printf("Ratio KY quarterly %f - Prod: %f and Capital: %f\n\n",KY*4, totalproduction, *totalassets);
    printf("Total labor = %f, Total Entrep labor = %f, Corp = %f, share_positive labor = %f\n\n",*totallabor, laborJ, *totallabor-laborJ,frac_employer/fracJ);
    printf("Share employers = %f %f %f \n",frac_employer_ability[0],frac_employer_ability[1],frac_employer_ability[2]);
    printf("MASS: U = %f, W = %f, E = %f\n", fracUi+fracUn, fracW, fracJ);
    printf("Fraction of E>nonE [6.7] %f\n",100*(fracJtoU + fracJtoWW)/fracJ);
    printf("Ratio median E/W = %f,  E/All = %f \n", mworthJ/mworthW, mworthJ/mworthall);
    printf("GINI(wealth) = %f - GINI(income) = %f, | std inc bus = %f, std inc work = %f, ratio = %f\n", GINIval, GINIearn, stdentinc, stdworkinc, stdentinc/stdworkinc);
    printf("Bankruptcy rate [0.75] %20.15f, rlv exit [10] %20.15f, rlv pop [0.164] %20.15f \n",100*(defaultrate/fracJ),100*((defaultrate/fracJ)/((fracJtoWW + fracJtoU)/fracJ)), 100*(defaultrate/(fracJ+fracUi+fracUn+fracW)));
    printf("U>E [2.2] %20.15f\n",100*(fracUntoJ + fracUitoJ)/(fracUi+fracUn));
    //printf("(UL: %20.15f, US: %20.15f) Unemploy. rate [6.4] %20.15f\n",fracUn, fracUi, 100*(fracUi + fracUn));
    printf("W>E 1 [1.07] %20.15f\n",transiWWtoE[0]/((fracWWtoJ)/fracW));
    printf("W>E 2 [0.89] %20.15f\n",transiWWtoE[1]/((fracWWtoJ)/fracW));
    printf("W>E 3 [1.07] %20.15f\n",transiWWtoE[2]/((fracWWtoJ)/fracW));
    printf("W>U 1 [1.07] %20.15f\n",transiWWtoU[0]/((fracWWtoUi)/fracW));
    printf("W>U 2 [0.89] %20.15f\n",transiWWtoU[1]/((fracWWtoUi)/fracW));
    printf("W>U 3 [1.07] %20.15f\n",transiWWtoU[2]/((fracWWtoUi)/fracW));
    printf("Neg income: %20.15f, Neg bus income: %20.15f \n",fracnegincome/fracJ, fracnegbusincome/fracJ);
    printf("Ratio AX/A: %f, Tot_assetJ: %f, CapitalJ/K: %f \n", ratioAXA, totalassetsJ/(*totalassets+totalassetsJ), capitalJ/(*totalassets+capitalJ));
    printf("Zero wealth fraction [7- 13] %20.15f\n",100*(distasset[0]));
    printf("flow(UtoJ) = %f, flow(UtononU) = %f, flow(JtoW) = %f, flow(WtononW) = %f\n",(fracUntoJ + fracUitoJ)/(fracUntoJ + fracUitoJ + fracWWtoJ), (fracUitoWW + fracUitoJ + fracUntoWW + fracUntoJ)/((fracUn + fracUi)), (fracJtoWW)/(fracJ), (fracWWtoJ + fracWWtoUi)/(fracW));
    printf("Fraction Ui/U [37] = %f\n", fracUi/(fracUi + fracUn));
    printf("IUR (Ui/W) [2.5] = %f\n", (fracUi/(fracW))*100);
    printf("-------TRANSITION Worker to entrepreneur / ABILITY--------\n");
    printf("\t very low \t middle \t very high\n");
    //for(tgrid = 0.0; tgrid < maxtheta; tgrid++)
    //{
    //    printf("%20.15f\t",aggWWtoJ[tgrid]/fracWWpability[tgrid]);
    //}
    //printf("\n");
    for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
    printf("%20.15f\t",transiWWtoE[tgrid]);
    }
    printf("\n");

    for(tgrid = 0.0; tgrid < maxtheta; tgrid++){
    printf("%20.15f\t",(necessityUnTHETA[tgrid]+necessityUiTHETA[tgrid]+necessityWWTHETA[tgrid])/fracJpability[tgrid]);
    }
    printf("\n");

    printf("--------- TRANSITION MATRIX BTW OCCUPATIONS ------\n");
    printf("  \t  U     \t W     \t E \n");
    printf("U \t %3.7f \t %3.7f \t %3.7f \n", (1- (fracUitoJ + fracUitoWW + fracUntoWW + fracUntoJ)/(fracUi+fracUn)), (fracUitoWW + fracUntoWW)/(fracUi+fracUn), (fracUntoJ + fracUitoJ)/(fracUi+fracUn));
    printf("W \t %3.7f \t %3.7f \t %3.7f \n", (fracWWtoUi)/fracW, (1 - (fracWWtoJ + fracWWtoUi)/fracW), (fracWWtoJ)/fracW);
    printf("E \t %3.7f \t %3.7f \t %3.7f \n", (fracJtoU)/fracJ, fracJtoWW/fracJ, (1-(fracJtoWW + fracJtoU)/fracJ));
    printf("---------------\n");
#endif




/*********************************/
/** PROJECTED METHOD SAVE FILES **/
/*********************************/

/******** SAVE IN A FILE FOR PROJECTED METHODS ******/
// (see top of the EUP.h to see how) //
#if SEA == 0
    saveDISTforPM(distJE, distJNE, distWWNE, distWWE, distUnE, distUnNE, distUiNE, distUiE);
#endif

#if SEA == 1
    saveDISTforPM(distJE, distJNE, distJEi, distJNEi, distWWNE, distWWE, distUnE, distUnNE, distUiNE, distUiE);
#endif





/****************************/
/** RETURN LIST OF MOMENTS **/
/****************************/

/**********  SET GENERATED MOMENT **********/
// global
generatedmoment[0] = KY;                    // Capital - output ratio
//if(generatedmoment[0] <= 2.65*1.05 && generatedmoment[0] >= 2.65*0.99){generatedmoment[0] = 2.65;} // acceptable bound

generatedmoment[1] = (defaultrate/fracJ)*100;   // default rate
if(generatedmoment[1] <= 0.57*1.05 && generatedmoment[1] >= 0.57*0.95){generatedmoment[1] = 0.57;} // acceptable bound

// labor market
generatedmoment[2] = 100*((fracJtoWW + fracJtoU)/(fracJ));       // Entrepreneur's exit rate
if(generatedmoment[2] <= 6.0*1.02 && generatedmoment[2] >= 6.0*0.98){generatedmoment[2] = 6.0;} // acceptable bound
generatedmoment[3] = 100*fracJ;       // Fraction of Entrepreneurs
generatedmoment[4] = 100*(fracUi + fracUn);       // Unemployment rate
generatedmoment[5] = 100*((fracUntoJ + fracUitoJ)/(fracUntoJ + fracUitoJ + fracWWtoJ));         // frac E previously U
if(generatedmoment[5] <= 19.4*1.025 && generatedmoment[5] >= 19.4*0.975){generatedmoment[5] = 19.4;}    // acceptable bound

// wealth income distribution
generatedmoment[6]  = (mworthJ/mworthW);        // ratio net worth
generatedmoment[7]  = (fracnegbusincome/fracJ)*10;        // zero net worth (*10) (because less precision)
if(generatedmoment[7] <= 1*1.08 && generatedmoment[7] >= 1*0.92){generatedmoment[7] = 1;} // acceptable bound

// transition btw occupations
generatedmoment[8]  = transiWWtoE[0]/((fracWWtoJ)/fracW);
generatedmoment[9]  = transiWWtoE[1]/((fracWWtoJ)/fracW);
generatedmoment[10] = transiWWtoE[2]/((fracWWtoJ)/fracW);
            





// PRINT STATISTICS FOR PAPER
#if indexPM == 0
    FILE *simulfile;

    char buffer_SS_statistics[150]; // The filename buffer.
    snprintf(buffer_SS_statistics, sizeof(buffer_SS_statistics), "OUTPUT/SS_statistics_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);

    simulfile = fopen(buffer_SS_statistics, "w");
        fprintf(simulfile,"fracE\tK_E\tK_Corp\tFracUtoE\tFracWtoE\tExitE\tBankruptcy\tUrate\tNecessity_share\tCorp_jobs\tE_jobs\tY\ttax_rate\tcost_over_gdp\twelfare\tC_welfaretot\tfracE1\tfracE2\tfracE3\tUtoE1\tUtoE2\tUtoE3\tUtoW1\tUtoW2\tUtoW3\tfracUtoW\tfracUi\tfracEi\tnbJ1\tnbW1\tnbU1\tnbJ2\tnbW2\tnbU2\tnbJ3\tnbW3\tnbU3\ttotal_labor\n");
        fprintf(simulfile,"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f",fracJ,capitalJ,totalcapital-capitalJ,(fracUntoJ + fracUitoJ)/(fracUi+fracUn),(fracWWtoJ)/fracW,(fracJtoWW + fracJtoU)/fracJ,100*(defaultrate/fracJ),100*(fracUi + fracUn),(necessityUn + necessityUi + necessityWW)/(fracUitoJ + fracUntoJ + fracWWtoJ),fracW,laborJ,totalproduction,*optimaltax,costSEA/totalproduction,welfaretot,C_welfaretot,fracJpability[0],fracJpability[1],fracJpability[2],aggUtoJ[0]/fracUpability[0],aggUtoJ[1]/fracUpability[1],aggUtoJ[2]/fracUpability[2],aggUtoWW[0]/fracUpability[0],aggUtoWW[1]/fracUpability[1],aggUtoWW[2]/fracUpability[2],(fracUntoWW + fracUitoWW)/(fracUi+fracUn),fracUi,fracJi,fracJpability[0], fracWWpability[0],fracUpability[0],fracJpability[1], fracWWpability[1],fracUpability[1],fracJpability[2],fracWWpability[2],fracUpability[2],*totallabor);
    fclose(simulfile);
#endif



/*************************************/
/** COMPUTE STATISTICS FOR 10 YEARS **/
/*************************************/



// The goal of this part is to simulate the economic dynamics of different group of new entrepreneurs.

// The first will apply this to three groups:
// group1: new entrepreneurs from paid employment,
// group2: from unemployment,
// group3: from necessity.
// for the moment group4 is useless.
#if SIMUL10_baseline == 1
    if(critprice < 0.01){
        
        printf("STARTING SIMULATION OF %d YEARS to recover survival rate and other statistics...\n", maxtimesimul);
        
        double *distnewJE_W, *distnewJNE_W, *distnewJE_U, *distnewJNE_U, *distnewJE_NEC, *distnewJNE_NEC;
        distnewJE_W         = (double *) calloc((ifulldimE), sizeof(double));
        distnewJNE_W        = (double *) calloc((ifulldimE), sizeof(double));
        distnewJE_U         = (double *) calloc((ifulldimE), sizeof(double));
        distnewJNE_U        = (double *) calloc((ifulldimE), sizeof(double));
        distnewJE_NEC       = (double *) calloc((ifulldimE), sizeof(double));
        distnewJNE_NEC      = (double *) calloc((ifulldimE), sizeof(double));


        // DEFINE THE GROUPS HERE //
        for(i = 0; i<maxgrid; i++) {
            for(t = 0; t<maxtheta; t++) {
                for(z = 0; z<maxfirmtype; z++) {
                    for(y = 0; y<maxtheta; y++) {
                        distnewJE_W[inxE(i,t,z)]    += mzinv[z]*(distWWEtoJE[inxW(i,t,y)]);
                        distnewJNE_W[inxE(i,t,z)]   += mzinv[z]*(distWWEtoJNE[inxW(i,t,y)]+distWWNEtoJNE[inxW(i,t,y)]);
                        distnewJE_NEC[inxE(i,t,z)]  += mzinv[z]*(necessityWWE[inxW(i,t,y)]);
                        distnewJNE_NEC[inxE(i,t,z)] += mzinv[z]*(necessityWWNE[inxW(i,t,y)]);
                    }
                    distnewJE_NEC[inxE(i,t,z)]  += mzinv[z]*(necessityUiE[inx(i,t)] + necessityUnE[inx(i,t)]);
                    distnewJNE_NEC[inxE(i,t,z)] += mzinv[z]*(necessityUiNE[inx(i,t)] + necessityUnNE[inx(i,t)]);
                    distnewJE_U[inxE(i,t,z)]     = mzinv[z]*(distUiEtoJE[inx(i,t)]);
                    distnewJNE_U[inxE(i,t,z)]    = mzinv[z]*(distUiEtoJNE[inx(i,t)] + distUiNEtoJNE[inx(i,t)]);
                }
            }
        }
        
    #if NOPOLICY == 1
        SIMULATION_ENT(distnewJNE_W, distnewJE_W, distnewJNE_U, distnewJE_U, distnewJNE_NEC, distnewJE_NEC, valueJNE, valueJE, valueUiE, valueUnE, valueWWE, valueWWNE, valueUnNE, valueUiNE, defaultJNE, collateralJNE, sloanJNE, saveresJNE, searchJNE_W, isaveJNE, collateralJE, saveresJE, searchJE_W, isaveJE, totalcapital, corporate_labor, 0);
    #endif
        
    #if SEA == 1
        SIMULATION_ENT(distnewJNE_W, distnewJE_W, distnewJNE_U, distnewJE_U, distnewJNE_NEC, distnewJE_NEC, valueJNE, valueJE, valueUiE, valueUnE, valueWWE, valueWWNE, valueUnNE, valueUiNE, defaultJNE, collateralJNE, sloanJNE, saveresJNE, searchJNE_W, isaveJNE, collateralJE, saveresJE, searchJE_W, isaveJE, defaultJNEi,collateralJNEi, sloanJNEi, saveresJNEi, searchJNEi_W, isaveJNEi, collateralJEi, saveresJEi, searchJEi_W, isaveJEi, valueJNEi, valueJEi, rateJNEi, totalcapital, corporate_labor, 0);
    #endif
        

        free(distnewJE_W);free(distnewJNE_W);free(distnewJE_U);free(distnewJNE_U);free(distnewJE_NEC);free(distnewJNE_NEC);
    }
#endif





// Group 1: new entrepreneurs coming from Unemployment, with downside risk insurance
// Group 2: new entrepreneurs coming from Unemployment, with subsidies
// Group 3: new entrepreneurs coming from Unemployment, without downside risk insurance neither subsidy.
// Group 4: counterfactual, new entrepreneur coming from Unemployment, with DRI, but same as in 3
// Group 5: counterfactual, new entrepreneur coming from Unemployment, with subsidy, but same as in 3
// Group 6: new entrepreneurs selected with the policy (residual between reform and baseline SS).

// Proceed as follows:
/*
1. save the counterfactual fraction, new entrepreneur coming from unemployment, without DRI neither subsidy.
2. for this fraction, give them the experiments, compare how their behavior differ from group 3.
3. evolution of group 1 and 2 are very different, because you have "composition effect" <-> not the same guys.
*/

// --> depending on the experiment, their is only three groups to be considered.

// Allow to compare //

#if writeSIMUL10 == 1
    FILE *writeUiEtoJE_bench, *writeUitoJNE_bench;
    // here we write the fraction of new entrepreneur coming from insured unemployment in the baseline model (Group. 3) //
    /**** UNEMPLOYED ****/
    writeUiEtoJE_bench=fopen(UiEtoJE_bench, "w"); setbuf ( writeUiEtoJE_bench, NULL );
    writeUitoJNE_bench=fopen(UitoJNE_bench, "w"); setbuf ( writeUitoJNE_bench, NULL );
    for(igrid=0;igrid<maxgrid;igrid++){
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
            for(int z = 0; z < maxfirmtype ; z++){
                fprintf(writeUiEtoJE_bench,"%20.15f\n",mzinv[z]*distUiEtoJE[inx(igrid,tgrid)]);
                fprintf(writeUitoJNE_bench,"%20.15f\n",mzinv[z]*(distUiEtoJNE[inx(igrid,tgrid)] + distUiNEtoJNE[inx(igrid,tgrid)]));
            }
        }
    }
    fclose(writeUiEtoJE_bench);
    fclose(writeUitoJNE_bench);
#endif


#if SIMUL10 == 1

//group1 = if(valueJNEi[] > valueXX[i] && valueXX[i] > valueJNE[i]) // <-- selected new entrepreneur with the insurance policy
//group1 = if(valueJNEs[] > valueXX[i] && valueXX[i] > valueJNE[i]) // <-- selected new entrepreneur with subsidy
//group2 = total INSURED      // i.e. total comming from Ui
//group2 = total subsidied    // i.e. total comming from Ui
//group3 = total newENT       // i.e. tout ceux qui come.
//for the moment group4 is useless.

double *distnewJE_Ui, *distnewJNE_Ui, *distnewJE_tot, *distnewJNE_tot, *distnewJE_baseline, *distnewJNE_baseline, *distnewJE_selected, *distnewJNE_selected;


distnewJE_Ui        = (double *) calloc((ifulldimE), sizeof(double));
distnewJNE_Ui       = (double *) calloc((ifulldimE), sizeof(double));
distnewJE_tot       = (double *) calloc((ifulldimE), sizeof(double));
distnewJNE_tot      = (double *) calloc((ifulldimE), sizeof(double));
distnewJE_baseline  = (double *) calloc((ifulldimE), sizeof(double));
distnewJNE_baseline = (double *) calloc((ifulldimE), sizeof(double));
distnewJE_selected  = (double *) calloc((ifulldimE), sizeof(double));
distnewJNE_selected = (double *) calloc((ifulldimE), sizeof(double));

if(critprice < 0.01){
    printf("STARTING SIMULATION OF %d YEARS...\n", maxtimesimul);

    // THE SAME FRACTION WHO ENTER ENTREPRENEUR AS UNEMPLOYED INSURED IN THE BENCHMARK ECONOMY [group 4 or 5] //
    readinput(distnewJE_baseline, ifulldimE, "OUTPUT/bench/distUiEtoJE.out");
    readinput(distnewJNE_baseline, ifulldimE, "OUTPUT/bench/distUitoJNE.out");


    // DEFINE THE GROUPS HERE //
    for(i = 0; i<maxgrid; i++) {
        for(t = 0; t<maxtheta; t++) {
            for(z = 0; z<maxfirmtype; z++) {

                // THOSE WHO ARE INSURED UNEMPLOYED, WHO BENEFITS FROM THE PROGRAM (GROUP 1 or 2 [group 3 if no policy is activated]) //
                distnewJE_Ui[inxE(i,t,z)] = mzinv[z]*distUiEtoJE[inx(i,t)];
                distnewJNE_Ui[inxE(i,t,z)] = mzinv[z]*(distUiEtoJNE[inx(i,t)] + distUiNEtoJNE[inx(i,t)]);
                // IF SEA == 1, then you will get that JEi and JNEi are always resulting from Ui :: so diff is in SIMUL_10y //

                // TOTAL (just to see) //
                distnewJE_tot[inxE(i,t,z)] = mzinv[z]*(distUnEtoJE[inx(i,t)]);
                distnewJNE_tot[inxE(i,t,z)] = mzinv[z]*(distUnNEtoJNE[inx(i,t)]  + distUnEtoJNE[inx(i,t)]);
                for(y = 0; y<maxyprod; y++) {
                    distnewJE_tot[inxE(i,t,z)] += mzinv[z]*(distWWEtoJE[inxW(i,t,y)]);
                    distnewJNE_tot[inxE(i,t,z)] += mzinv[z]*(distWWNEtoJNE[inxW(i,t,y)] + distWWEtoJNE[inxW(i,t,y)]);
                }
            }
        }
    }


    // GROUP WHO IS SELECTED (group 6) // (just need to change either baseline or Ui in the next function by selected.
    for(i = 0; i<maxgrid; i++) {
        for(t = 0; t<maxtheta; t++) {
            for(z = 0; z<maxfirmtype; z++) {
                distnewJE_selected[inxE(i,t,z)] = max((distnewJE_Ui[inxE(i,t,z)] - distnewJE_baseline[inxE(i,t,z)]),0);
                distnewJNE_selected[inxE(i,t,z)] = max((distnewJNE_Ui[inxE(i,t,z)] - distnewJNE_baseline[inxE(i,t,z)]),0);

//                if(distnewJNE_selected[inxE(i,t,z)] > 0.000000001) {
//                    printf("ERROR JNE %20.15f %20.15f %20.15f", distnewJNE_selected[inxE(i,t,z)],distnewJNE_baseline[inxE(i,t,z)],distnewJNE_Ui[inxE(i,t,z)]);getchar();
//                }
//                if(distnewJE_selected[inxE(i,t,z)] > 0.000000001) {
//                    printf("ERROR JE %20.15f %20.15f %20.15f", distnewJE_selected[inxE(i,t,z)],distnewJE_baseline[inxE(i,t,z)],distnewJE_Ui[inxE(i,t,z)]);getchar();
//                }
            }
        }
    }


    #if NOPOLICY == 1
    SIMULATION_ENT(distnewJNE_baseline, distnewJE_baseline, distnewJNE_selected, distnewJE_selected, distnewJNE_Ui, distnewJE_Ui, valueJNE, valueJE, valueUiE, valueUnE, valueWWE, valueWWNE, valueUnNE, valueUiNE, defaultJNE, collateralJNE, sloanJNE, saveresJNE, searchJNE_W, isaveJNE, collateralJE, saveresJE, searchJE_W, isaveJE, totalcapital, corporate_labor,1);
    #endif

    #if SEA == 1
    SIMULATION_ENT(distnewJNE_baseline, distnewJE_baseline, distnewJNE_selected, distnewJE_selected, distnewJNE_Ui, distnewJE_Ui, valueJNE, valueJE, valueUiE, valueUnE, valueWWE, valueWWNE, valueUnNE, valueUiNE, defaultJNE, collateralJNE, sloanJNE, saveresJNE, searchJNE_W, isaveJNE, collateralJE, saveresJE, searchJE_W, isaveJE, defaultJNEi,collateralJNEi, sloanJNEi, saveresJNEi, searchJNEi_W, isaveJNEi, collateralJEi, saveresJEi, searchJEi_W, isaveJEi, valueJNEi, valueJEi, rateJNEi, totalcapital, corporate_labor,1);
    #endif

}
#endif





/*********************************/
/** SAVE IN FILES BY OCCUPATION **/
/*********************************/

#if logprint == 1

/**** ENTREPRENEURS ****/
distoutJ=fopen(distfileJ, "w"); setbuf ( distoutJ, NULL );
fprintf(distoutJ,"%5s\t%5s\t%20s\t%20s\t%20s\n","z", "THETA", "ASSET", "distJE", "distJNE");
for(int zgrid = 0; zgrid < maxfirmtype; zgrid++) {
    for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
        for(igrid=0;igrid<maxgrid;igrid++){
            fprintf(distoutJ,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\n",zgrid, tgrid, grid[igrid], distJE[inxE(igrid,tgrid,zgrid)],distJNE[inxE(igrid,tgrid,zgrid)]);
        }
    }
}
fclose(distoutJ);




/**** ENTREPRENEURS INCOME DIST ****/
distoutJ=fopen(distfileJ, "w"); setbuf ( distoutJ, NULL );
fprintf(distoutJ,"%5s\t%5s\t%20s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZETA", "ASSET", "z", "distJE", "busincJE", "incJE", "distJNE", "busincJNE", "incJNE", "sloanJNE", "distJNEstay", "distJEstay");
for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
for(zgrid=0;zgrid<maxfirmtype;zgrid++){
for(igrid=0;igrid<maxgrid;igrid++){
    for(z=0; z<maxfirmtype; z++) {
    
        if(sloanJNE[inxE(igrid, tgrid, zgrid)] > 0){indicatorAX = 1.0;}else{indicatorAX = 0.0;}
        
        fprintf(distoutJ,"%5d\t%5d\t%20.15f\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, grid[igrid], z, distJE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z],
        
        (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJE[inxE(igrid, tgrid, zgrid)])),
        
        (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJE[inxE(igrid, tgrid, zgrid)]) + rstar*(collateralJE[inxE(igrid, tgrid, zgrid)] - grid[igrid])),
        
        distJNE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z],
        
        (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJNE[inxE(igrid, tgrid, zgrid)]) - (rateJNE[inxE(igrid, tgrid, zgrid)]-1.0)*(sloanJNE[inxE(igrid, tgrid, zgrid)])),
        
        (1.0 - corptax)*(profitE(statez[z],tgrid,collateralJNE[inxE(igrid, tgrid, zgrid)]) - (rateJNE[inxE(igrid, tgrid, zgrid)]-1.0)*(sloanJNE[inxE(igrid, tgrid, zgrid)])) + rstar*(1.0 - indicatorAX)*(grid[igrid]-collateralJNE[inxE(igrid, tgrid, zgrid)]),
        
         sloanJNE[inxE(igrid, tgrid, zgrid)],
         
         mzprob[zgrid][z]*(distJNE[inxE(igrid,tgrid,zgrid)] - distJNEtoUnE[inxE(igrid,tgrid,zgrid)] - distJNEtoUiE[inxE(igrid,tgrid,zgrid)] - distJNEtoUiNE[inxE(igrid,tgrid,zgrid)] - distJNEtoUnNE[inxE(igrid,tgrid,zgrid)] - distJNEtoWWE[inxE(igrid,tgrid,zgrid)] - distJNEtoWWNE[inxE(igrid,tgrid,zgrid)]), mzprob[zgrid][z]*(distJE[inxE(igrid,tgrid,zgrid)] - distJEtoUnE[inxE(igrid,tgrid,zgrid)] - distJEtoWWE[inxE(igrid,tgrid,zgrid)] -  distJEtoUiE[inxE(igrid,tgrid,zgrid)] - distJEtoWWNE[inxE(igrid,tgrid,zgrid)]));
}}}}
fclose(distoutJ);

#if SEA == 1
distoutJ=fopen(distfileJi, "w"); setbuf ( distoutJ, NULL );
fprintf(distoutJ,"%5s\t%5s\t%20s\t%5s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZETA", "ASSET", "z", "distJEi", "incJEi", "distJNEi", "incJNEi");
for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
    for(zgrid=0;zgrid<maxfirmtype;zgrid++){
        for(igrid=0;igrid<maxgrid;igrid++){
        fprintf(distoutJ,"%5d\t%5d\t%20.15f\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, grid[igrid],z, distJE[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z],(1.0 - corptax)*(profitE(statez[z],tgrid,collateralJEi[inxE(igrid, tgrid, zgrid)]) - rstar*(collateralJEi[inxE(igrid, tgrid, zgrid)]- grid[igrid])), distJNEi[inxE(igrid,tgrid,zgrid)]*mzprob[zgrid][z],(1.0 - corptax)*(profitE(statez[z],tgrid,collateralJNEi[inxE(igrid, tgrid, zgrid)]) - rstar*(collateralJNEi[inxE(igrid, tgrid, zgrid)] - grid[igrid])));
}}}
fclose(distoutJ);
#endif


/**** UNEMPLOYED ****/
distoutU=fopen(distfileU, "w"); setbuf ( distoutU, NULL );
fprintf(distoutU,"%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ASSET", "INCOMEUn", "distUnE", "distUnNE",  "INCOMEUi", "distUiE", "distUiNE");
for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
for(igrid=0;igrid<maxgrid;igrid++){
fprintf(distoutU,"%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid], minpar, distUnE[inx(igrid,tgrid)],distUnNE[inx(igrid,tgrid)], prod[tgrid]*rhostar*wstar*(1.0 - taxrate), distUiE[inx(igrid,tgrid)],distUiNE[inx(igrid,tgrid)]);
}}
fclose(distoutU);



/**** WORKER ****/
distoutWW=fopen(distfileW, "w"); setbuf ( distoutWW, NULL );
fprintf(distoutWW,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\n","THETA", "PROD", "ASSET", "INCOME", "distWWE", "distWWNE");
for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
for(int ygrid = 0; ygrid < maxyprod; ygrid++){
for(igrid=0;igrid<maxgrid;igrid++){
    fprintf(distoutWW,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, ygrid, grid[igrid], prod[tgrid]*yprod[ygrid]*wstar*(1.0 - taxrate)+rstar*grid[igrid],distWWE[inxW(igrid,tgrid, ygrid)],distWWNE[inxW(igrid,tgrid, ygrid)]);
}}}
fclose(distoutWW);


/**** WEALTH ****/
distoutwealth=fopen(distfilewealth, "w");setbuf ( distoutwealth , NULL );
fprintf(distoutwealth,"%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","ASSET", "all", "J","WW","Ui","Un");
for(igrid=0;igrid<maxgrid;igrid++){
fprintf(distoutwealth,"%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",grid[igrid],distasset[igrid],distassetJ[igrid],distassetW[igrid],distassetUi[igrid],distassetUn[igrid]);
}
fclose(distoutwealth);



/**** FIRM SIZE ****/
distoutfirm=fopen(distfilesize, "w"); setbuf ( distoutfirm , NULL );
fprintf(distoutfirm, "%5s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\n","TYPE","PROD","ASSET","investJNE","distJNE","investJE","distJE");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(igrid=0;igrid<maxgrid;igrid++){
for(zgrid=0; zgrid<maxfirmtype; zgrid++){
    fprintf(distoutfirm,"%5d\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",zgrid, tgrid, igrid, (sloanJNE[inxE(igrid, tgrid, zgrid)] + collateralJNE[inxE(igrid, tgrid, zgrid)]), distJNE[inxE(igrid, tgrid, zgrid)], collateralJE[inxE(igrid, tgrid, zgrid)], distJE[inxE(igrid, tgrid, zgrid)]);
}}}
fclose(distoutfirm);



// PRINT TRANSITION DISTRIBUTION J TO WW AND WW TO J
distoutfile=fopen(flowJtoWW, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"%5s\t%5s\t%20s\t%20s\t%20s\n","TYPE", "PROD", "ASSET", "JEtoWW", "JNEtoWW");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(igrid=0;igrid<maxgrid;igrid++){
for(zgrid=0;zgrid<maxfirmtype;zgrid++){
    fprintf(distoutfile,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\n",zgrid, tgrid, grid[igrid], distJEtoWWE[inxE(igrid, tgrid, zgrid)] + distJEtoWWNE[inxE(igrid, tgrid, zgrid)], distJNEtoWWE[inxE(igrid, tgrid, zgrid)] + distJNEtoWWNE[inxE(igrid, tgrid, zgrid)]);
}}}
fclose(distoutfile);

distoutfile=fopen(flowWWtoJ, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"%5s\t%5s\t%20s\t%20s\n","PROD", "ASSET", "WWtoJ", "JtoWW");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(ygrid=0; ygrid<maxyprod; ygrid++){
for(igrid=0;igrid<maxgrid;igrid++){
    fprintf(distoutfile,"%5d\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid], distWWEtoJE[inxW(igrid, tgrid, ygrid)] + distWWEtoJNE[inxW(igrid, tgrid, ygrid)], distWWNEtoJNE[inxW(igrid, tgrid, ygrid)]);
}}}
fclose(distoutfile);

distoutfile=fopen(flowUiEtoJ, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"%5s\t%5s\t%20s\t%20s\n","PROD", "ASSET", "UiEtoJ", "UiNEtoJ");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(igrid=0;igrid<maxgrid;igrid++){
fprintf(distoutfile,"%5d\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid], distUiEtoJE[inx(igrid, tgrid)] + distUiEtoJNE[inx(igrid, tgrid)], distUiNEtoJNE[inx(igrid, tgrid)]);
}}
fclose(distoutfile);


// PRINT NECESSITY SHARE DISTRIBUTION //
// for Unemployed
distoutfile=fopen(distnecessity, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"%5s\t%5s\t%20s\t%20s\n","PROD", "ASSET", "Un", "Ui");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(igrid=0;igrid<maxgrid;igrid++){
fprintf(distoutfile,"%5d\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid], necessityUnE[inx(igrid, tgrid)] + necessityUnNE[inx(igrid, tgrid)], necessityUiE[inx(igrid, tgrid)] + necessityUiNE[inx(igrid, tgrid)]);
}}

fclose(distoutfile);
// for Worker
distoutfile=fopen(distnecessityW, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"%5s\t%5s\t%20s\t%20s\n","PROD", "ASSET", "WWE", "WWNE");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(ygrid=0; ygrid<maxyprod; ygrid++){
for(igrid=0;igrid<maxgrid;igrid++){
    fprintf(distoutfile,"%5d\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid], necessityWWE[inxW(igrid, tgrid, ygrid)], necessityWWNE[inxW(igrid, tgrid, ygrid)]);
}}}
fclose(distoutfile);


// To find threshold value at which enter entrepreneurship, conditional of receiving job and business opportunity //

#if computethreshold == 1
// 1. compute the expected value of being entrepreneur //
double *expvalueE, *threshold, stop;
expvalueE = (double *) calloc((ifulldim), sizeof(double));
threshold = (double *) calloc((maxtheta*maxyprod), sizeof(double));

for(tgrid = 0; tgrid<maxtheta; tgrid++){
for(igrid = 0; igrid < maxgrid; igrid++){
for(zgrid = 0; zgrid < maxfirmtype; zgrid++){
    expvalueE[inx(igrid, tgrid)] += mzinv[zgrid]*valueJNE[inxE(igrid, tgrid, zgrid)];
}}}

// 2. find the threshold value (of wealth) for any y, t //
for(tgrid = 0; tgrid<maxtheta; tgrid++){
    for(ygrid = 0; ygrid < maxyprod; ygrid++){
        igrid = 0;
        stop = 0;

        for(igrid = 0; igrid < maxgrid; igrid++){
            if(expvalueE[inx(igrid, tgrid)] >= valueWWNE[inxW(igrid, tgrid, ygrid)]){ // not this one, look at the UiNE
                threshold[ygrid*maxtheta+tgrid] = grid[igrid];
                break;
            }
        }
    }
}


distoutfile=fopen(thresholdfile, "w");
setbuf ( distoutfile , NULL );
fprintf(distoutfile,"%5s\t%5s\t%20s\n","THETA", "Y", "THRESHOLD");
for(tgrid=0;tgrid<maxtheta;tgrid++){
for(ygrid=0; ygrid<maxyprod; ygrid++){
fprintf(distoutfile,"%5d\t%5d\t%20.15f\n",tgrid, ygrid, threshold[ygrid*maxtheta+tgrid]);}}
fclose(distoutfile);
free(expvalueE);
free(threshold);
#endif

#endif






